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In the last decade, level-set methods have been actively developed for applications in 
image registration, segmentation, tracking, and reconstruction. However, the development 
of a wide variety of level-set PDEs and their numerical discretization schemes, coupled 
with hybrid combinations of PDE terms, stopping criteria, and reinitialization strategies, 
has created a software logistics problem. In the absence of an integrative design, current 
toolkits support only specific types of level-set implementations which restrict future 
algorithm development since extensions require significant code duplication and effort. 
In the new NIH/NLM Insight Toolkit (ITK) v4 architecture, we implemented a level-set 
software design that is flexible to different numerical (continuous, discrete, and sparse) 
and grid representations (point, mesh, and image-based). Given that a generic PDE is 
a summation of different terms, we used a set of linked containers to which level-set 
terms can be added or deleted at any point in the evolution process. This container-based 
approach allows the user to explore and customize terms in the level-set equation 
at compile-time in a flexible manner. The framework is optimized so that repeated 
computations of common intensity functions (e.g., gradient and Hessians) across multiple 
terms is eliminated. The framework further enables the evolution of multiple level-sets 
for multi-object segmentation and processing of large datasets. For doing so, we restrict 
level-set domains to subsets of the image domain and use multithreading strategies 
to process groups of subdomains or level-set functions. Users can also select from a 
variety of reinitialization policies and stopping criteria. Finally, we developed a visualization 
framework that shows the evolution of a level-set in real-time to help guide algorithm 
development and parameter optimization. We demonstrate the power of our new 
framework using confocal microscopy images of cells in a developing zebrafish embryo. 
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1. INTRODUCTION 

The automated identification of anatomical structures found in 
medical and microscopy images is an important step in any 
imaging-based quantitative analysis pipeline. Large variations in 
image quality arising from differences in acquisition protocols, 
anisotropic point-spread functions, and image noise complicate 
the task of automated image analysis tools. To overcome the dis- 
advantages associated with simple heuristic methods, a class of 
methods for contour evolution known as geometric active con- 
tour models, or level-sets, have been actively developed (Osher 
and Sethian, 1988; Sethian, 1999). Level-sets have become a 
preferred method for addressing a number of image science prob- 
lems (Tsai and Osher, 2004) including denoising (Rudin et al., 
1992; Dibos and Koepfler, 1998; Vese and Osher, 2004), reg- 
istration (Vemuri et al, 2000, 2003; Droske and Ring, 2006), 
segmentation (Caselles et al, 1995; Malladi et al, 1995; Leventon 
et al, 2000; Cremers et al, 2006), tracking (Dufour et al., 2005; 
Dydenko et al, 2006; Dzyubachyk et al., 2010), and surface recon- 
struction (Zhao et al, 2001; Nilsson et al., 2005). Level-set meth- 
ods represent the presumptive boundary C of an object of interest 
as the zero level-line of a higher dimensional implicit function 



C(t) = l(x, y)\<ty(x, y, t) = 0}, also called the level-set function. 
For example, the boundary C can be arbitrarily initialized along 
with an initial level-set function <$>o(x, y) constructed as a signed 
distance function to C. In Figure 1A, we plot the zero curve of a 
level-set function initialized with a square. Then, the evolution 
of the level-set function to match the true object boundary is 
governed by setting speed functions or via the minimization of 
energy functionals (Figures 1B-D). In a basic formulation, the 
evolution equation of the level-set function can be specified as 
follows: 

+P|Vch| = 0 (1) 

of 

$(x,y, 0) = rj>o(x, y) 

The function F is called the speed function and depends on 
the image data I as well as (j). An advantage of the level-set 
method, especially for medical imagery, is its natural ability to 
incorporate information on object shape, texture, and color dis- 
tribution into the segmentation process. Level-sets also avoid 
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FIGURE 1 | (A-D) Iterative stages of active contour (red) evolution using 
level-sets for the segmentation of objects (cells in this case) in images. 
The underlying image is a 2D confocal microscopy section showing 
nuclei (H2B:GFP) in the zebrafish inner ear. (E) The level-set function is 



typically defined over the entire 2D image domain. Different 
narrow-band representations are shown in (F) Whitaker, (G) Shi, and 
(H) Malcolm. The Whitaker method uses 5 layers, Shi uses 2, and 
Malcolm uses 1 layer. 



the problem of explicit parameterization of the object bound- 
ary, a problem with parametric active contour approaches such 
as snakes (Kass et al., 1988), thus providing flexibility in the 
segmentation of objects with topological changes, cusps, and cor- 
ners. Furthermore, the same level-set formulation applied on 2D 
images naturally extends to an N-dimensional image. Given the 
large number of level-set methods already developed, we refer 
to (Osher and Fedkiw, 2004) for a complete exposition of the 
level-set calculus. We next discuss the major classes of level-sets 
methods and how they impact the development of a generic 
software system. 

2. BACKGROUND 

In general, there are two main classes of the level-set methodology 
that arise from using a variational approach that minimizes an 
energy functional: (i) edge-based and (ii) region-based methods. 

2.1. EDGE-BASED LEVEL-SETS 

The goal of edge-based approaches is to evolve the contour until 
it finds an object edge where its speed is gradually reduced to 0. 
An edge-based level-set based on mean curvature motion is given 
by (Caselles et al, 1993): 



<ty(x, y, 0) 



g(V/)]V(|>| | div 

<t>o(*, y) 



(2) 



where v > 0 is a constant. The update equation is constructed 
from the image intensity (I), gradient (VI), and edge function 
g(V7). The edge function g(VJ) is designed to have values close to 



0 at the boundary and large values elsewhere. The zero level-curve 
of this level-set moves with a speed of g(I)\ V4>| (div + 
and stops on the object boundary where g approaches zero. The 
constant term ensures the level-set always advances toward the 

object boundary even when the curvature (div ^^jL-^) becomes 

negative, i.e., ^div^y^-^ + v ) remains positive. Another pop- 
ular example of an edge-based level-set uses the strong image 
gradient at the object boundary to slow down or stop the zero 
curve (Malladi et al., 1993, 1994): 



3(|> 

f t = IV+I 
<\>(x,y, 0) = $ 0 (x,y) 



v + 



(Ml - M 2 ) 



(|VG*/|-M 2 ) (3) 



where v is a constant, G is a Gaussian function, and Mi and 
M 2 are the maximum and minimum values of the image gradi- 
ent (|VG*J|). In these level-sets, parts of the contour reach the 
boundary and cross-over before the rest of the contour catches 
up. In order to prevent the contour from over-shooting the edge 
so that it remains trapped in a minimum along the boundary, 
Kichenassamy et al. (1995) and Caselles et al. (1995) indepen- 
dently proposed the geodesic active contour formulation. Here, 
the underlying energy is representative of the contour length in a 
Riemannian space with a metric induced by the image intensity: 



3f v V |V<)>| 



+ vg(|VJ|) 



(4) 
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= g(I)\V$\div + VgO) • V<|> + v#(|VJ|) 

<fy(x,y, 0) = 4>oO,y) 

The first (curvature) and last (propagation) terms in Equation 4 
approach zero at the object boundary. The middle (advection) 
term ensures that contour motion is always directed toward the 
boundary. 

Many additional variants of edge-based concepts have been 
published. An abstract representation common to all edge-based 
partial differential equation (PDE) is as follows: 



dt 



-aA(x) ■ V(|> - PP(x)| Vcp| + yZ(x)k:|V(j)| (5) 



<K*,y, 0) = <\>o(x,y) 

where A is an advection term, P is a propagation (expansion) 
term, and Z is a spatial modifier term for the mean curvature k. 
The scalar constants a, p\ and y weight the relative influence of 
each of the terms on the movement of the interface. Based on this 
prevalent model in the early 2000s, level-sets were implemented 
in the NIH/NLM Insight Toolkit (ITK) v3. As we show next, 
this model is inadequate in its representation of region-based 
level-sets. 

2.2. REGION-BASED LEVEL-SETS 

Region-based level-sets segment the image into objects based on 
region statistics (rather than just object edges) of intensity, tex- 
ture, or color values. For example, the region mean intensity for 
foreground (ci) and background {c^) is a popularly used statis- 
tic for defining an energy functional F whose minimization leads 
to an optimal segmentation of the foreground (Chan and Vese, 
1999, 2001): 

F(c l: c 2 , <t>) = / U(x) - ci) 2 dx 

JlnsideiQ 

+ / (J(x) - c 2 ) 2 dx + v ■ Area{C) (6) 

JOutside(C) 

+\x, ■ Length(C) 

Minimizing the energy functional using the standard gradient 
descent (or steepest descent) method leads to the following 
gradient-flow equation: 



dl 



i* ||||J -v-^cr-d) 2 (7) 



+ ^2(i-c 2 y 



Cl = 



f n I(x, y).H(${x, y))dxdy 
Jo. H($(x, y))dxdy 

f Q I(x,y).(l-H(<$>{x,y)))dxdy 
f Q (l-H(<\>(x,y)))dxdy 



$(x,y, 0) = <$>o(x,y) 



where H is the Heaviside function, S e is the smoothed Dirac-Delta 
function, and u, and v regularize the curve length and fore- 
ground area, respectively. To account for image inhomogeneities 
and large gradients that may be present in an image, piecewise- 
smooth extensions were proposed in Tsai et al. (2001); Vese and 
Chan (2002); Li et al. (2008). In these extensions, local intensity- 
functions are defined in place of using constants c\ and c 2 and the 
energy functional was defined in terms of the fit with respect to 
these functions. Nevertheless, the final equation form is similar to 
Equation 7. 

Within the last decade, region-based techniques, as well as 
graph-partitioning-based active contours (Sumengen et al, 2004; 
Sumengen and Manjunath, 2006), have emerged. These new 
methods do not ascribe to the same generic PDE (Equation 5) 
used in ITKv3 framework. The addition of region-based meth- 
ods in the ITKV3 framework required us to duplicate significant 
amounts of code (Mosaliganti et al., 2009a,b,c). Although edge 
and region-based methodologies arise from different strategies, 
our new software design and implementation in the ITK v4 com- 
bines the two hierarchies seamlessly with minimal code duplica- 
tion. New terms can also be easily added and require no changes 
in the evolution and update classes. 

2.3. NARROW-BAND METHODS 

One of the primary disadvantages with using level-set methods 
for image segmentation is that they are slow and memory- 
intensive. The level-set function is typically discretized on the 
entire image grid to hold floating-point values (Figure IE and 
Movie_Dense_Sl.mov) although only the position of the zero 
level-curve is of primary interest. Therefore, a key develop- 
ment in reducing computational cost has been the emer- 
gence of narrow-band algorithms (Whitaker, 1998; Malcolm 
et al, 2008; Shi and Karl, 2008). These methods evolve the 
level-set function in a layer around the zero level-set (sparse- 
field) alone as opposed to its solution on the entire image 
domain (dense). While the Whitaker method (Whitaker, 1998) 
(Figure IF and Movie_Whitaker_S2.mov) uses 5 layers around 
the zero level-curve, the Shi and Karl (2008) (Figure 1G and 
Movie_Shi_S3.mov) and Malcolm et al. (2008) (Figure 1H 
and Movie_Malcolm_S4.mov) implementations use 2 and 1, 
respectively. Together, they provide a significantly faster and 
less memory-intensive implementation although the trade-off 
is that the solution is only maintained around the zero level- 
curve. They may also lead to a different local solution closer 
to the initialization contour compared to the dense case. In 
current software systems, narrow-band implementations derive 
from the dense case and thus use an image-based repre- 
sentation and do not fully exploit the speed and memory 
enhancement possible. In our new software design, we derive 
the sparse level-set formulations in a separate hierarchy that 
represents the level-set function as linearized lists of pix- 
els with floating-point values. The pixels map onto a label 
image stored in run-length encoded format. By introducing 
suitable narrow-band base classes, we have allowed for future 
enhancements and new narrow-band representations to be intro- 
duced without much effort or code-duplication in the ITKv4 
framework. 
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2.4. REINITIALIZATION AND STOPPING CRITERION 

During the evolution of the level-set function, the differential 
movement of different level-curves can cause the function to 
develop gradients that are too steep and/or too flat. This neg- 
atively impacts the numerical stability in subsequent iterations 
leading to unpredictable outcomes. Therefore, one needs to per- 
form a distance reinitialization step at the end of every few itera- 
tions, such that the zero level-curve location remains unchanged. 
One way to perform such an initialization is to evolve the PDE to 
a steady-state according to Sussman et al. (1994): 

9(b 

+ s^«(cp 0 )(|Vct>| - 1) = 0 (8) 

of 

<$>(x,y, 0) = §o(x,y) 

where 4>o represents the level-set function before the reinitial- 
ization. The end result will be a signed distance function to the 
interface (4>o = 0). Another approach is to solve the Eikonal 
equation using fast-marching methods (Peng et al, 1999): 

IV+I = 1 (9) 

<K*, y, t) = 0 <^> <j)o(x, y) = 0 

In this method, the signed distance function is used to fix the 
level-set function in a narrow band around the zero curve as 
boundary conditions. Fast marching methods are then used to 
solve the Eikonal equation. 

In another development, Li et al. defined a new term in the 
equation that automatically maintains the level-set function to a 
signed distance function (Li et al., 2005). 

P(4>) = ( ^(|V<H - ifdxdy (10) 

The addition of this energy followed by the steepest descent 
flow leads to the addition of a new term ( ArJ> — div fj^fj) = 

div ^(1 — )V(f>j. In regions where |V(()| > 1, this term causes 
an outward diffusion of level-curves. Similarly, where |V((>| < 1, 
there is an inward diffusion of level-curves. The term approaches 
zero when | V(j)| = 1 leading to a solution of the Eikonal equation. 

Similar to the reinitialization problem, the stopping criterion 
also presents several choices to the user. The evolution of a level- 
set function is typically halted by a threshold set on the number of 
iterations (N), and/or by assessing the reduction in the variational 
energy, and/or by assessing the change in the level-set function, 
and/or by checking to see if the level-set has reached certain 
pre-set boundary points. In order to let users explore the above 
strategies and possibly develop new ones, we implemented reini- 
tialization and stopping base classes that serve as plugins into the 
level-set evolution framework in ITK v4. This provides complete 
flexibility to the user in developing a customized implementation 
of level-sets without restrictions from the design. 

2.5. MULTI-OBJECT AND MULTIPHASE METHODS 

In biomedical and biological image analysis applications, one 
is often interested in segmenting more than a single object (of 



the same or different kind) from a given image. This espe- 
cially happens when the objects to be segmented are adjacent 
to each other and the delineation of one object automatically 
affects the neighboring object. For example, in MRI imagery, 
white and gray matter regions in the brain have a common 
interface. Moreover, spatial constraints often specify relation- 
ships between classes of level-sets. For example, microscopy 
images depict thousands of cells in terms of their plasma mem- 
branes, nuclear membranes, organelles, and proteins. Here, the 
nucleus and organelles are always contained within the mem- 
brane. In such situations, concurrently interacting and evolving 
level-sets provide optimal and efficient image segmentation solu- 
tions. Multi-object methods (Dufour et al., 2005; Mosaliganti 
et al., 2009a,c) use a unique level-set function per anatomical 
structure in the image. Overlaps among the level-set functions 
can be penalized or encouraged through a special term placed 
in each of the level-set update equations. Another strategy is to 
use multiphase methods, where different phases of the set of 
level-set functions (+ and — parts) encode a unique anatomi- 
cal structure (Vese and Chan, 2002). While multiphase strategies 
work well for a small number of objects, multi-object strategies 
are more suitable for the case of a larger number of objects. 
In our new software framework, we enable the solution for a 
large system of level-set equations. No restriction is placed on 
the type of level-sets (dense or sparse). Each level-set function 
may evolve according to a different update equation specified 
by the user. Furthermore, the number of level-set functions in 
the system may dynamically change during iterations (say, from 
cell division or from a cell entering or exiting the field-of-view 
during tracking). The individual terms of the update equations 
can also dynamically change. Thus, accounting for such flex- 
ibility in the design allows one to fully exploit the level-set 
techniques. 

3. METHODS 

The goal of the new inheritance framework is to permit the 
implementation of a variety of level-set methods with minimal 
code-duplication and allow end-users to explore and customize 
the method to their data. The new framework is designed to 
be modular, and therefore, reduces the effort needed to expand 
the level-set code base. An overview of the level-set algorithm 
with components of our framework is provided in Figure 2. 
Corresponding class inheritance diagrams for the components are 
shown in Figures 3, 5. We elaborate on the components of the 
new framework next. 

3.1. DOMAIN REPRESENTATION 

While the level-set function has typically been implemented 
on underlying image grids, it also theoretically extends to 
unstructured grids or surface meshes and continuous repre- 
sentations. In the popular ITK v3 framework, for example, 
an image representation was used throughout the code which 
impedes the development of mesh-based level-set segmenta- 
tion methods. Additionally, in a continuous representation, the 
level-set function is not discretized on a grid but instead rep- 
resented analytically in terms of linear combination of base 
functions [e.g., radial basis functions (Bernard et al., 2007; 
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Algorithm overview 

Initialize level-sets 
Compute domain mapping 
Begin iteration 
a) Initialize term parameters 
Compute update 
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FIGURE 2 | Components of the level-set framework. (A) An overview of 
the basic level-set algorithm is provided. (B) Individual level-set domains are 
assumed to be subsets of the image domain. (C) The domain is partitioned 
into a label-map that specifies the interacting level-set ids and the 



corresponding subregion. (D) Level-set functions, equations, and, terms are 
stored in a container format. During iteration, the computation of the level-set 
updates are multithreaded either in terms of processing individual domains or 
individual level-set functions. 



Gelas et al., 2007), splines (Bernard et al., 2009)]. Such an 
underlying representation may still use the same energy for- 
mulation and equation terms. Thus, modularizing the level- 
set code base such that the term and equation classes do 
not directly depend on the domain representation is essen- 
tial for building a common framework. We next detail how 
the level-set function is defined to account for all such 
representations. 

3.2. LEVEL-SET FUNCTION 

Our code is heavily templated to allow for a great deal of 
flexibility in level-set methods in an efficient manner. In 
our new framework, we define an abstract level-set func- 
tion base class (itk : : LevelSetBase) inheriting from 
itk : : DataObj ect based on 4 template parameters 
(Figure 3A). 

template< class TInput, unsigned int 
VDimension, typename TOutput, 

class TDomain > class LevelSetBase : 
public DataObj ect 

where TInput defines the input type where the level-set func- 
tion will be evaluated, VDimension is dimension of the input 
space, TOutput is the returned type when evaluating the level-set 
function (for the general case when it is not a scalar), TDomain 
is the discretization of the level-set function (e.g., ImageBase or 
QuadEdgeMesh). In the case of an image representation, this spe- 
cializes as only the first three parameters and inherits from the 
itk: : ImageBase class: 



template< class TInput, unsigned int 
VDimension, class TOutput > class 
LevelSetlmage : public LevelSetBase< 
TInput, VDimension, TOutput, ImageBase< 
VDimension > > 

While this definition accounts for an image representation 
(continuous or discrete), we further specialize into a discrete 
image representation with TInput as itk: : Index to enable 
the traditional image-discretized implementation of the level-set 
method. 

template< typename TOutput, unsigned int 
VDimension > 

class DiscreteLevelSetlmage : public 
LevelSetImage< Index< VDimension >, 
VDimension, TOutput > 

All of the above level-set function classes implement specific 
member functions for returning the level-set value [<\>(x, y)], gra- 
dient (V((>), Hessian ( V 2 ())), Laplacian (rj)^ + (|>w)> gradient norm 
(|V(|>|), and mean curvature (k = div(j^j-)) given its underly- 
ing representation (continuous or discrete image or mesh). Thus, 
the level-set equation, term, and evolution classes are indepen- 
dent of the underlying domain representation which facilitates 
the implementation of a wide variety of level-set methods. 

The level-set quantities are stored as specific instances of a 
templated class itk :: LevelSetBase :: DataType<T>. 
The template parameter T determines whether the stored 
member variable is a scalar, vector, or tensor. The class 
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itlclevelSetDomainPartition 
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itlcLevelSetDomainPartition 
Base< TImage > 



itk::LevelSetDomainPartition 
Image< TImage > 



itlclevelSetDomainPartition 
ImageWithKdTree< TImage > 



itlcDataObject 



X 



itk::Object 



itk::LevelSetBase< 
TInput, VDimension, 
TOutput TDomain > 



itlc:LevelSetImage 
< Index< VDimension 
>, VDimension, TOutput > 
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rerm 
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itk::LevelSetEqua 
Base< TInput, TL 


tionTerm 

evelSetContainer > 



itlc:DiscreteLevelSetImage 
TOutput VDimension > 



itlcLevelSetDenselmage 
< TImage > 



itlcLevelSetSparselmage 
< TOutput, VDimension > 
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itlcShiSparseLevelSetlmage 
< VDimension > 



FIGURE 3 | Inheritance diagrams for (A) Domain partitioning classes, (B) Level-set function, and (C) Level-set equation term classes in the ITK v4 
framework. 



contains three member variables that store the quan- 
tity name, value, and a Boolean flag indicating if it had 
already been computed. The level-set quantities are then 
collectively instantiated inside a special structure called 
itk: :LevelSetBase: : LevelSetDataType. When 
computing the update in a level-set equation (Figure 2A), the 
different terms in the equation repetitively compute the level-set 
quantities. Therefore, by passing the LevelSetDataType object 
among various terms, we ensure that level-set quantities are 
computed exactly once, cached, and reused. This significantly 
reduces the computations involved for calculating each term and 
the update overall. 

The discrete image representation itk : : Discre 
teLevelSetlmage is then specialized into the dense 
(itk : : LevelSetDenselmage) and sparse (itk:: 
LevelSetSparselmage) cases, which in turn was spe- 
cialized into three sparse representations (itk: :Whi taker 
SparseLevelSetlmage, itk: : ShiSparseLevelSet 
Image, and itk: :MalcolmSparseLevelSetImage) 
(Figures 1E-H, 3B). The sparse representation (also some- 
times called as narrow-band) discretizes the level-set function 
only around the zero level-curve. The number of layers 
and their update schemes differ among the three variants. 
A layer is defined as a map data-structure that links image 
indices to level-set function values. Each layer is associated 
with a unique id and the set of all the layers are stored as 
another map data- structure linking the id values to the layer. 
Additionally, all the sparse-field classes store a label-map 
(run-length image encoding) of the layers which significantly 
reduces the amount of memory used over a regular image 
(Lehmann, 2007). 



3.2.1. Image to level-set adaptors 

For the convenience of end-users, adaptors for converting from 
binary images to level-set function objects were developed in ITK 
v4. The supplied binary image could be the output of a prepro- 
cessing segmentation pipeline. We consider the dense level-set 
function as well as the three types of sparse-field representations: 
(1) Whitaker, (2) Shi, and (3) Malcolm. The adaptor base class 
is templated over the input image type and the output level-set 
function type. It contains member variables for the input image 
and the computed level-set function representation and a pure 
virtual member function for computing the level-set function. 

template< class TInputlmage, class 
TLevelSet > class 

Binary ImageToLevelSetlmageAdap tor Base : 
public Object 

The derived class itk: : Binary ImageToLevelSetlmage 
Adaptor uses the partial template specialization mechanism 
for defining the virtual member function corresponding to each 
level-set representations. 

3.3. RESTRICTED LEVEL-SET DOMAINS 

The large computational costs associated with level-set methods is 
a particularly severe problem in multi-level set implementations. 
Typically, depending on the computing environment, one cannot 
evolve more than tens of level-set functions without running into 
memory problems. To address this problem, our new framework 
handles level-set functions defined on a subset of the input image 
domain (Q,). The interaction among level-set functions is limited 
to those functions whose domains overlap (Figures 2B,C). A 
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helper base class (LevelSetDomainPartitionBase) is 
used to define the location and size of the level-set domains 
relative to f2 (Figure 3A). The domain of each level-set function 
is stored in a vector (m_LevelSetRegionVector). The class 
contains two pure virtual methods AllocateListDomain ( ) 
and PopulateListDomain ( ) that allocate and pop- 
ulate a new data-structure (mesh or image) depending 
on the underlying grid. Each grid point stores a list of 
the active level-set function ids. The class is specialized to 
LevelSetDomainPartitionlmage in the case of a dis- 
crete image grid and to LevelSetDomainPartitionMesh 
for the case of meshes. For the case when there are thousands of 
level-sets, populating a list image by checking overlap at each pixel 
is time-consuming. Therefore, we further specialized into a class 
itk : : LevelSetDomainPartitionlmageWithKdTree. 
This class uses a iCd-tree data structure that contains the cen- 
troids of the level-set domains. The Kd- tree is used to query 
nearby level-set functions at each pixel and check for overlap. 
This enables the simultaneous evolution of thousands of level-set 
functions thereby expanding the applicability of level-set proce- 
dures to tracking large numbers of objects and in large images. 
Note that there is an initial overhead associated with building the 
Kd-tree that can be avoided for cases involving a small number 
of level-set functions. 

The output of the helper class is to define an image of list 
pixel types. Each list holds the ids of level-set functions active 
at that pixel. By using another helper class (itk: : Level 
SetDomainMapImageFilter), the list image is further clus- 
tered into sub-regions based on pixel types. Each sub-region has 
the same set of active level-set ids and can therefore be iter- 
ated upon for computing the update efficiently (Figure 2C). The 
itk : : LevelSetDomainMapImageFilter is a member of 
the itk: : LevelSetContainerBase. 

As an illustrative example.Figure 4 shows the application of 
the domain partitioning technique on a sample 2D confocal 
image of the zebrafish ear. In Figure 4A, we show closely-packed 
nuclei with varying intensities and noise. Figure 4B shows rect- 
angular domains placed around each nuclei that define the extent 
of the initialized level-set function. The rectangular domains 
are subdivided to encode a unique level-set interaction within. 
Figure 4C shows the final level-set function and Figure 4D shows 
an overlay of the segmentation on the raw image. The small 
rectangular domain of each level-set does not influence its evolu- 
tion and final segmentation but dramatically reduces the memory 
utilized and total running-time. 

3.4. TERMS 

The level-set equation is a weighted sum of terms. Each term 
is a simple function of the level-set (e.g., gradient, Hessian, 
divergence, etc.), the input image (intensity, gradient, etc.), or 
involving the other level-set functions concurrently evolving in 
the system. As mentioned before, level-set computations occur 
repeatedly in each term and therefore values can be computed, 
cached, and stored in the level-set object and shared among 
the various terms. The term base class implements functions 
[Evaluate ( . ) ] for computing the contribution from a 
term toward the level-set update. It also contains member 



functions for initializing/updating term parameters before 
[initialize ( ) ]) and after [Update ( ) ] an iteration. While 
these functions are used for updating parameters in a dense 
level-set method, concurrent updating [UpdatePixel ( . ) ] 
is carried out in a sparse level-set method. This is because only 
small incremental changes occur in the level-zero curve position 
in a sparse level-set method. Concurrent updating ensures that 
term parameters stay current and cuts down on computations 
in between iterations. This further enhances the speedup from 
using sparse level-set methods. The computed update is finally 
scaled by a weight data member (m_Coef f icient). The 
term class also stores a level-set function identifier to supply the 
computed update to the evolution of the corresponding level-set 
function in the container. The term base class is templated using 
two parameters, namely, the input image type and a container of 
level-set functions (Figure 3C): 

template< class TInputlmage , class 
TLevelSetContainer > 

class LevelSetEquationTermBase : public 
Obj ect 

Different types of terms arising from edge-based and 
region-based level-set methods such as the propaga- 
tion, Laplacian, advection, curvature, and region-based 
terms described in Equations 4 and 7 derive directly from 
LevelSetEquationTermBase: 

template< class TInput, class 

TLevelSetContainer> 

class LevelSetEquationLaplacianTerm : 
public LevelSetEquationTermBase< TInput, 

TLevelSetContainer > 

3.5. CONTAINER-BASED DESIGN 

Consider a concurrently evolving system of level-set functions 
{4>i , 4>2 , • • • , 4>m} with input image I, and domain £2. The generic 
level-set equation (l/j) consists of a weighted summation of a 
number of different terms Ty. 

[7i: ^M = g ai .. ri . (UW M i) 

TT 3<t>2(/>) t a \ x. \M \ 

U 2 : — — = 2^ a 2j • T 2j (I, {(t),}™ ,) 

T i=i 



U M ■ dr = 2^ a M; • T Mj(l {<)><}, = i) 

The set of terms that constitute each equation may vary across all 
the M equations. Additionally, the number of concurrently evolv- 
ing level-set functions in the system may dynamically change. 
To enable such flexibility, we used containers to store level-set 
function objects, equation objects, and their constitutive terms 
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FIGURE 4 1 Cell segmentation using multiple level-sets defined on 
subsets of the image domain. (A) 2D confocal image of cell 
nuclei in the zebrafish ear. (B) Square domains defined around 
approximate centers of each cell. (C) Visualization of the 



(Figure 2D). In each case, the container holds elements along 
with an id that is common across the level-set function, equation, 
and term containers. They are implemented as class templates, 
which allows a great flexibility in the types supported as elements. 
The container manages the storage space for its elements and pro- 
vides member functions to access them, either directly or through 
iterators (reference objects with similar properties to pointers) or 
by using the id. 

3.5.1. Level-set container 

In the case of multiple level-set functions, a special container 
class (itk: : LevelSetContainerBase) was templated 
over two input parameters, namely, identifier type of the level-set 
function and the level function type (dense or sparse). The 
class provides iterators (const and non-const) and member 
functions for adding/removing individual level-set objects 
[AddLevelSet ( ) and GetLevelSet ( ) ] by using the 
level-set identifier. Member functions also allow for copying 
and performing logical comparisons among container objects. 
A domain map filter (m_DomainMapFilter) is also instanti- 
ated that describes how the domain is split among the different 
level-set objects (see section 3.3). 

template< class TIdentifier, class 
TLevelSet > 

class LevelSetContainerBase : public 
Obj ect 

Using the partial template specialization mechanism, the class 
itk: : LevelSetContainer is specialized depending on 
the type of the level-set function. The derived class imple- 
ments a member function for allocating new memory and 
copying individual level-set functions [Copylnf ormation 
AndAllocate ( . ) ]. For example, the following class special- 
izes the definition for the dense level-set method: 

template< class TIdentifier, class TImage > 
class LevelSetContainer< TIdentifier, 
LevelSetDenseImage< TImage > > : 




overlapping level-set functions. (D) Final segmentation obtained using 
the geodesic active contour method. The level-set used advection, 
propagation, curvature, and overlap penalty terms with weights of 
3, 3, 1, and 1000, respectively. 



public LevelSetContainerBase< TIdentifier, 
LevelSetDenseImage< TImage > > 

3.5.2. Term containers 

For each level-set function defined by its identifier 
(CurrentLevelSetld), a term container class is instan- 
tiated that holds all the terms in its level-set equation. Using 
AddTerm ( ) and PushTerm ( ) member functions, new 
terms are introduced/removed from the level-set equation. 
While calculating the update to a level-set, the term container 
object is iterated and each term is evaluated and summed 
up to yield the update. Thus, member functions for initial- 
izing [initializeParameters ( ) ], updating individual 
term parameters [Update ()], and evaluating all terms 
[Evaluate ( ) ] are provided. 

template< class TInputlmage, class 
TLevelSetContainer > 

class LevelSetEquationTermContainer : 
public Object 

The class is templated over the input image type and the level-set 
container type. The class is not derived further since there is no 
further specialization of the individual terms. In a given term con- 
tainer, the end-user can combine edge-based and region-based 
terms for customizing the level-set method which is yet another 
aspect of flexibility enabled by our framework. 

3.5.3. Level-set equation container 

For a system of level-set equations, each equation ([/,) cor- 
responds to a unique term container. The class LevelSet 
EquationContainer is templated over the type of the 
term container. LevelSetEquationContainer holds the individual 
term containers for all the level-set functions. A unique term 
container can be referenced by the level-set function identifier (i) 
through GetEquation ( ) member function. New equations 
can be added to the container using AddEquation ( ) . While 
calculating the update to a level-set, the equation container object 
is iterated and each term container is in turn updated Update 
InternalEquationTerms ( ) . Thus, member functions 
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< TEquationContainer. 
LevelSetDenseImage< 
TImage > > 



itlclevelSetEvolutionBase 
< TEquationContainer, LevelSetDense 
Image< TImage > > 



itlclevelSetEvolution 
< TEquationContainer, 
LevelSetDenseImage< TImage > 



itlc:Object 



itlclevelSetEvolutionBase 

< TEquationContainer, TLevelSet > 



< TEquationContainer, 
MalcolmSparseLevelSetlmage 

< VDimension > > 

< TEquationContainer, 
WhttakerSparseLevelSetlmage 

< TOutput, VDimension > > 

< TEquationContainer, 
ShiSparseLevelSetlmage 

< VDimension > > 



itlclevelSetEvolutionBase 
< TEquationContainer, MakolmSparse 
LevelSetImage< VDimension > > 



itlclevelSetEvolution 
< TEquationContainer, 
MalcolmSparseLevelSetlmage 
< VDimension > > 



itlclevelSetEvolutionBase 

< TEquationContainer, WhitakerSparse 

LevelSetImage< TOutput, VDimension > 



itlclevelSetEvolution 
< TEquationContainer, 
WhitakerSparseLevelSetlmage 
< TOutput, VDimension > > 



itlclevelSetEvolutionBase 
< TEquationContainer, ShiSparse 
LevelSetImage< VDimension > > 


< 


itlclevelSetEvolution 
< TEquationContainer, 
ShiSparseLevelSetlmage 
< VDimension > > 









FIGURE 5 | Inheritance mechanism for level-set evolution classes. 



for initializing (initializeParameters ( ) ) and updating 
individual term parameters (Update ( ) ) are provided. 

template< class TTermContainer > 

class LevelSetEquationContainer : public 

Object 

3.6. LEVEL-SET EVOLUTION 

The class itk : : LevelSetEvolutionBase implements the 
iterations in the level-set method in the Evolve ( ) member 
function. At the beginning of each iteration, update buffers 
are allocated [AllocateUpdateBuf f er ( ) ] and term 
parameters are initialized [initializelteration ( ) ]. 
Assuming the stopping criterion function is not satisfied, the 
iteration proceeds by computing the updates in the equation 
container, which in turn calls the individual term contain- 
ers. The updates are computed by stepping through smaller 
domains where the active level-set function ids are known (from 
section 3.3). This avoids the costly step of checking at each 
pixel if a particular level-set is active. Based on the updates 
computed, the overall time-step for evolving the level-set is 
determined [ComputeTimeStepForNextlteration ( ) ]. 
Based on the time-step, the level-set functions are updated 
[UpdateLevelSets ( ) ] and the term parameters are then 
updated [UpdateEquations ( ) ]. Based on the change in the 
level-set functions, the global change is accumulated and set as 
input to the stopping criterion function. This function (described 
next) determines whether the next iteration continues. The 
base class is templated over the equation container type and the 
level-set function type. 

template< class TEquationContainer, class 
TLevelSet > 

class LevelSetEvolutionBase : public 
Object 

Dense and sparse-strategies are different in terms of how the 
update buffers are setup. Sparse methods visit only the narrow 



band layers while the dense method visits every pixel in the 
level-set domain. Therefore, as before we use partial template 
specialization to specialize the implementation depending on the 
level-set representation. 

Figure 6 shows the application of our framework to the auto- 
mated segmentation of cells from 3D confocal images of a 
developing zebrafish embryo expressing fluorescent proteins. The 
images show nuclei in magenta and membranes in green. The 
high-resolution images are of dimensions 1024 x 1024 x 72 
with a pixel sampling of 0.2 \im x 0.2 u,m x 1.0 |im. Using the 
Chan and Vese (1999, 2001) region-based terms with weights set 
to 1 (see section 2.2 and Equation 7) and an overlap penalty 
term of 1000, we evolved a total of 946 level-sets (485 cells and 
461 nuclei) (Dufour et al., 2005). Level-sets were initialized as 
small spheres with a radius of 3.0 |xm based on cell centers 
that were previous identified by using a Hough transform for 
identifying spherical objects. Figures 6A,B show the XY and XZ 
planes. Membrane and nuclei segmentations are contoured in yel- 
low and red, respectively. Figure 6C shows the domain partitions 
from the nuclei level-sets alone. Figure 6D shows the a three- 
dimensional rendering of the segmentation after 100 iterations of 
the system. 

The concurrent segmentation of all cells and correspond- 
ing cell nuclei in ITKv3 is intractable given that each level-set 
function would occupy the same domain as the image and 
946 level-set functions would not fit in memory. Secondly, 
the level-sets would have to be segmented independently one 
function at a time. This would lead to a running time of 
many hours. Thirdly, level-set functions cannot interact with 
each other. Figure 6E highlights these differences between ITKv3 
and the new framework in regards to the segmentation prob- 
lem. In order to make meaningful comparisons, we performed 
domain partitioning in ITKv3 apriori by using region-of-interest 
filters outside the level-set API. While a single level-set is 
evolved at a time, several intermediate processing filters such 
as those involved in computing level-set reinitialization are 
intrinsically multithreaded. Thus, ITKv3 execution, is partially 
multithreaded. 
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FIGURE 6 | Demonstration of cell segmentation in large 3D 
nuclear+membrane confocal microscopy data of a developing zebrafish 
ear at 24 h post-fertilization. The images depict cell membranes in green 
(membrane-targetted mCitrine) and nuclei in magenta (H2B-tdTomato). The 
images are of dimensions 1024 x 1024 x 72 with a pixel sampling of 0.2 |im 
x 0.2 |im x 1.0 p.m. A total of 485 whole cells and 461 nuclei were 
automatically identified by seeding level-sets at centers first identified by a 
Hough transform. (A) XY and (B) XZ sections showing membrane and nuclei 



segmentations in yellow and red contours, respectively. (C) A 3D 
randomly-colored rendering of the individual level-set domains. (D) A 3D 
rendering showing the segmentations all the ear cells with the tri-planar 
cross-sections of the image data. (E) Performance comparison between the 
old ITKv3 and the new ITKv4 frameworks, the number of threads used (T), 
with domain-discretization (DD), and accounting for level-set overlaps. The 
columns indicate if the framework can process multiple level-sets, if they 
allow interaction, and their running times, respectively. NA, not available 



In ITKv4, domain partitioning happens automatically and the 
user exactly specifies the number of threads to use. We provided 
running times for the procedure (~21 min). DD refers to domain 
partitioning and T refers to the number of threads used. We also 
show running times without overlap (~19min) and with the 
ITKv3 framework (~41 min). Thus, our method improves signif- 
icantly over ITKv3 in speed of computation even after external 
domain discretization as well as ease of use and available options. 
Finally, we also evolve geometrically interacting level-sets that 
register a small increase in running time. 

3.7. STOPPING CRITERION 

The iterations in the level-set method end when the 
functions finish evolving (converge) or when they satisfy 



some user-defined criterion. A specialized stopping criterion 
class (itk : : StoppingCriterionBase) is implemented 
and checked at the end of each iteration [isSatis 
fied()]. Typically, the method is terminated when there 
is no appreciable change in the level-set function or in 
the variational energy being minimized (itk : : Level Se 
tEvolutionStoppingCriterion). In some applications, 
it is more useful to set a higher-limit on the iteration num- 
ber (itk: : LevelSetEvolutionNumberOf Iterations 
StoppingCriterion). In other applications, the arrival at 
a certain fraction of user-defined boundary points might be 
more meaningful. Thus, the new design enables many such 
approaches by allowing users to choose the appropriate stopping 
criteria. 
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3.8. USER-INTERACTION 

In order to enable communication among filters and with graph- 
ical user-interfaces (GUIs), we used the itk: :EventObject 
classes. These provide a standard coding interface for sending 
and receiving messages indicating the initiation/termination of 
processes and modification of filters. At the end of each itera- 
tion of the level-set method, an event [iterationEvent ( ) ] 
is triggered. This allows users to use observers (callbacks) on 
the itk: : LevelSetEvolution object to execute a com- 
mand [AddObserver ( const EventObject & event, 
Command * ) ] . Examples of downstream operations may 
include updating a progress bar in a visualization pipeline or in 
an interactive level-set method where the user can seed or termi- 
nate level-sets. This user-interaction implementation follows the 
subject/observer design pattern. 

3.9. VISUALIZATION 

To enable the real-time visualization and recording of level- 
set function evolution, we provide a pipeline linking our 
level-set framework to The Visualization Toolkit (VTK). The 
process of visualization consists of first mapping a level-set 
function to a VTK dataset. The base class for converting 
a itk : : LevelSetFunction object is derived from 
itk : : ProcessObj ect and templated over the level-set 
function type. 

template< class TLevelSet > 

class LevelSetToVTKImageDataBase : public 

ProcessObj ect 

Given the dense and three sparse-field representations of level- 
set functions implemented (Figures 1E-H), we specialize a class 
itk: : LevelSetToVTKImageData for each type. In the 
dense level-set, the underlying level-set function is sampled on 
the image grid and passed to a VTK image defined on the same 
grid. In the sparse cases, the label-maps storing the locations of 
the sparse layers are converted into an ITK image and then to 
a VTK image. Having mapped the level-set function to a VTK 
image, the second step consists of developing suitable visualiza- 
tion of the level-set function using an image actor, renderer, and 
a rendering window that is annotated with the current iteration 
number. For a dense level-set function, one way is to visualize 
the level-set function with a simple colormap (itk: :Visual 
izelmageLevelSet). Another option, only available in a 
two-dimensional case, is to visualize the level-set function as a ele- 
vation map by using a surface mesh representation (itk: : VTKV 
isualize2DLevelSetAsElevationMap) (Figure IE). At 
each point (x, y) in the grid, the level-set function is eval- 
uated [z = <$>(x, y)] to build a set of interconnected vertices. 
The mesh is color-mapped based on the height (z). A third 
option is to visualize the level-curves (user-defined isoval- 
ues) of the function (itk: : VTKVisualizelmageLevel 
SetlsoValues) (Figure ID). For this, the marching squares 
(in 2D), marching cubes (in 3D) algorithm are used to 
extract contours/meshes from the image and then color-mapped. 
For the sparse representations, visualization consists of ren- 
dering the layers around the zero curve of the function 



using a colormap (itk: : vtkVisualize2DSparse-Level 
SetLayersBase). Each layer is mapped to a unique user- 
defined color (Figures 1F-H). Like before, the class is specialized 
using partial template specialization for the three sparse rep- 
resentations because they contain different numbers of layers. 
Evolution of level-sets using the dense (Supplementary Movie 
SI) and the three narrow-band representations (Supplementary 
Movies S2-S4) can be visualized in real-time using the inter- 
action classes described in Section 3.8 and the visualization 
classes described above. The level-sets were initialized as a square 
and segment the cell nuclei image shown in Figure 1 for 100 
iterations. 

4. IMPLEMENTATION AND DISCUSSION 

We chose to implement our framework using C++ in the 
NIH/NLM Insight Toolkit (ITK) v4 because of its solid software 
engineering practices, permissive license, community-support, 
and its popularity in the biomedical imaging community. The 
level-set classes, examples, and associated tests have been inte- 
grated into the ITKv4 toolkit. Detailed instructions for down- 
loading ITK are available at: http://www.itk.org/Wiki/ITK/Git. 
Instructions for compiling and installing ITK on all common 
computer systems are available at: http://www.itk.org/ITK/help/ 
documentation.html. 

After compiling and installing ITK, users can navi- 
gate to the downloaded source directory and find a total 
of 123 level-set classes at the following location: ITK/ 
Modules/Segmentation/LevelSetsv4/include. The level- 
set visualization classes that work with the Visualization 
Toolkit (VTK) are located at: ITK/Modules/Segmentation/ 
LevelSetsv4Visualization/include. Documentation for the 
level-set classes detailing the member functions, variables, 
and inheritance hierarchy is available online at: http://www. 

itk.org/Doxygen44/html/group ITKLevelSetsv4.html Each of 

the level-set classes is subject to unit tests individually as well as 
integrative tests. A total of 67 tests using simple datasets have 
been created to test the framework automatically. The tests are 
located at: ITK/Modules/Segmentation/LevelSetsv4/test. The 
tests are automatically compiled when building ITK and also 
provide further examples of code usage and API to new users. 
Users can run the tests from the compiled binary directory. The 
results of the automatic testing are reported on hundreds of 
computers across the world on a daily basis at: http://open.cdash. 
org/index.php?project=Insight. The results indicate whether 
the tests successfully configured, compiled, and executed on the 
remote system and also provide running times. 

Additional examples of the level-set API used to generate 
the results for this article (Figure 1 and Supplementary Movies 
S1-S4) have been provided online in the ITKExamples repository: 
https://github. com/InsightSoftwareConsortium/ITKExamples. 
The repository includes 2D test image data and information for 
compiling and running our code. Other examples of ITK APIs 
and level-set code are also available online in a third repository 
at:http://www.itk.org/Wiki/ITK/Examples#ImageSegmentation. 
Finally, as open-source developers, we provide continuous 
feedback and troubleshoot problems reported by users on the 
ITK mailing lists:http://www.itk.org/ITK/help/mailing.html. 
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Thus, our framework addresses a long-standing challenge in 
the image analysis community for developing a state-of-the art 
ITKv4 segmentation and analysis toolkit. Our framework achieves 
the following goals: 

4.1. SCALABILITY TO LARGE DATASETS AND EXTENSIBILITY TO NEW 
METHODS 

Current medical and microscopic imaging modalities feature 
rich datasets that requires the continuous development of new 
algorithms and performance optimization strategies to handle 
large data. In our framework, we combined the two major 
classes of edge-based and region-based methods seamlessly which 
should promote the development and usage of hybrid meth- 
ods. We refactored sparse-field implementations to enable the 
addition of Whitaker, Malcolm, and Shi methods for improving 
performance. We enabled the implementation of multi level- 
set methods to handle multiple objects in images and also 
implemented domain partitioning techniques for improved per- 
formance. Thus, our framework supports large datasets, incorpo- 
rates all the major level-set technologies, and enables the research 
and development of new ideas in the near future. 

4.2. ADAPTABILITY TO USER-CUSTOMIZATIONS AND DATASET 
VARIATIONS 

Our implementation allows the user to customize all aspects of 
the level-set method. The user can optimize the number of level- 
sets to use and specify how they interact. Users can specify the 
exact terms in the level-set update equations and choose the reini- 
tialization and stopping criteria to use after each iteration. The 
large number of customization options available means that the 
user is better equipped to handle the varying challenges presented 
by different biomedical datasets. Users can also seed or terminate 
level-set functions, or manipulate the terms and term weights at 
the beginning of each iteration in the system. This paves the way 
for the exploration and development of real-time level-set systems 
in the future. 

4.3. EFFICIENT PERFORMANCE BY MULTITHREADING, SPARSE-FIELD, 
AND DOMAIN PARTITIONING STRATEGIES 

The three sparse-field strategies significantly improve level-set 
performance on large datasets. Domain partitioning significantly 
reduces the memory utilized per level-set especially when the 
objects occupy a small portion of the large input space. By 
caching and reusing level-set quantities in the term calcula- 
tions (see Section 3.2), we further improve performance. We 
also implemented multithreaded execution of the iterations to 
take advantage of multicore processor systems. Multithreading 
is accomplished by allocating threads to iterate on (i) different 
regions in the domain partitioned image or (ii) across individual 
level-sets. This choice is decided at run-time depending on the 
maximum number of threads available to the system, the num- 
ber of level-sets initialized by the user, and the number of domain 
partitions present. 

5. FUTURE WORK 

Currently, there is no other comprehensive framework or toolkit 
similar to the ITK framework in terms of implementing a 



multitude of level-set technologies in a scalable and efficient 
manner. However, upon further testing we do find that ITKv3 per- 
forms 2-2.5 x faster for very basic segmentation using just one 
thread. This speed difference appears to result from the use of 
image iterators in ITKv3 compared with GetPixelQ and SetPixelQ 
methods in ITKv4 so future work could see a speed improve- 
ment by transitioning to image iterators. In the future, we plan to 
develop user- friendly APIs for developers familiar with ITKv3 and 
GUIs that will benefit new users. For the purpose of cell segmen- 
tation in microscopy, we have already begun developing a image 
analysis GUI called GoFigure2 (www.gofigure2.org). Currently, 
level-set implementations and applications have been restricted 
to discretized image domains. However, the level-set calculus is 
quite generic and can be extended to unstructured grids like 
meshes and analytically to continuous representations. Mesh seg- 
mentations are useful, for example, to segment surfaces. In the 
microscopy domain, this could serve to identify gene expression 
patterns on fly embryo surfaces or morphogen gradients on the 
tissue surfaces. Continuous representations can potentially give 
rise to faster level-sets on high-resolution image datasets. Real- 
time segmentation and tracking, and human-computer interac- 
tive applications are another major problem domain that could 
benefit from level-set appproaches. Real-time manipulation of 
level-sets in terms of the number of functions, their individual 
terms, and interactions can now be accomplished by our adaptive 
design. 
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